home *** CD-ROM | disk | FTP | other *** search
-
- { PROGRAM PI1.PAS
- SEPT 2, 1985
- ( IMPROVED VERSION OF PI.PAS WHICH WAS SLOW AND HAD SOME ERRORS )
-
- THIS PROGRAM COMPUTES THE DIGITS OF PI USING THE ARCTANGENT FORMULA
- (1) PI/4 = 4 ARCTAN 1/5 - ARCTAN 1/239
- IN CONJUNCTION WITH THE GREGORY SERIES
-
- (2) ARCTAN X = SUM (-1)^N*(2N + 1)^-1*X^(2N+1) N=0 TO INFINITY.
-
- SUBSTITUTING INTO (2) A FEW VALUES OF N AND NESTING WE HAVE,
-
- PI/4 = 1/5[4/1 + 1/25[-4/3 + 1/25[4/5 + 1/25[-4/7 + ...].].]
-
- - 1/239[1/1 + 1/239^2[-1/3 + 1/239^2[1/5 + 1/239^2[-1/7 +...].].]
-
- USING THE LONG DIVISION ALGORITHM, THIS ( NESTED ) INFINITE SERIES CAN BE
- USED TO CALCULATE PI TO A LARGE NUMBER OF DECIMAL PLACES IN A REASONABLE
- AMOUNT OF TIME. A TIME FUNCTION IS INCLUDED TO SHOW HOW SLOW THINGS
- GET WHEN N IS LARGE. IMPROVEMENTS CAN BE MADE BY CHANGING THE SIZE OF
- THE ARRAY ELEMENTS HOWEVER IT GETS A BIT TRICKY.
-
- A LITTLE HISTORY
- ----------------
-
- IN AUGUST, 1949, PROFESSOR JOHN VON NEUMANN USED THIS FORMULA TO
- CALCULATE PI TO OVER 2000 DECIMAL PLACES ON THE ENIAC COMPUTER.
- THE CALCULATION WAS COMPLETED OVER THE LABOR DAY WEEKEND WITH THE
- COMBINED EFFORTS OF FOUR ENIAC STAFF MEMBERS EACH WORKING EIGHT-HOUR
- SHIFTS TO ENSURE CONTINUOUS OPERATION OF THE ENIAC.
-
- SOME YEARS AGO I REQUESTED INFORMATION ON PI FROM THE ENCYCLOPEDIA
- BRITANNICA RESEARCH SERVICE. I RECEIVED A REPORT GIVING THE ABOVE
- HISTORICAL ACCOUNT PLUS A LISTING OF THE 2000 DIGITS.
-
- IT WAS THIS LISTING THAT ENABLED ME TO CHECK THE PROGRAM AND KEEP
- MY SANITY.
-
- HAVE FUN
- CINO HILLIARD
- [72756,672] }
- { PROGRAM PI1.PAS }
- Type
- TimeString = string[10];
- function time: TimeString;
- type
- regpack = record
- ax,bx,cx,dx,bp,di,si,ds,es,flags: integer;
- end;
- var
- recpack: regpack;
- ah,al,ch,cl,dh,dl: byte;
- hour,min,sec,hds: string[$2];
- BEGIN
- ah := 2;
- with recpack do
- begin
- ax := ah shl 8 + al;
- end;
-
- intr($1A,recpack);
- with recpack do
- begin
- str(cx shr 8,hour);
- str(cx mod 256,min);
-
- str(dx shr 8,sec);
- str(dx mod 256,hds);
- end;
- time := hour+':'+min+':'+sec+':'+hds;
- END;
-
- VAR B,C,V,P1,S,K,N,I,J,Q,M,M1,X,R,D:INTEGER;
- P,A,T:ARRAY[0..5000] OF INTEGER;TI:STRING[20];
- LABEL 10;
- CONST F1=5;
- CONST F2=239;
- PROCEDURE DIVIDE(D:INTEGER);
- BEGIN
- R:=0;
- FOR J:=0 TO M DO
- BEGIN
- V:= R*10+P[J];
- Q:=V DIV D;
- R:=V MOD D;
- P[J]:=Q;
- END;
- END;
- PROCEDURE DIVIDEA(D:INTEGER);
- BEGIN
- R:=0;
- FOR J:=0 TO M DO
- BEGIN
- V:= R*10+A[J];
- Q:=V DIV D;
- R:=V MOD D;
- A[J]:=Q;
- END;
- END;
- PROCEDURE SUBT;
- BEGIN
- B:=0;
- FOR J:=M DOWNTO 0 DO
- IF T[J]>=A[J] THEN T[J]:=T[J]-A[J] ELSE
- BEGIN
- T[J]:=10+T[J]-A[J];
- T[J-1]:=T[J-1]-1;
- END;
- FOR J:=0 TO M DO
- A[J]:=T[J];
- END;
- PROCEDURE SUBA;
- BEGIN
- FOR J:=M DOWNTO 0 DO
- IF P[J]>=A[J] THEN P[J]:=P[J]-A[J] ELSE
- BEGIN
- P[J]:=10+P[J]-A[J];
- P[J-1]:=P[J-1]-1;
- END;
- FOR J:= M DOWNTO 0 DO
- A[J]:=P[J];
- END;
- PROCEDURE CLEARP;
- BEGIN
- FOR J:=0 TO M DO
- P[J]:=0;
- END;
- PROCEDURE ADJUST;
- BEGIN
- P[0]:=3;
- P[M]:=10;
- FOR J:=1 TO M-1 DO
- P[J]:=9;
- END;
-
- PROCEDURE ADJUST2;
- BEGIN
- P[0]:=0;
- P[M]:=10;
- FOR J:=1 TO M-1 DO
- P[J]:=9;
- END;
-
- PROCEDURE MULT4;
- BEGIN
- C:=0;
- FOR J:=M DOWNTO 0 DO
- BEGIN
- P1:=4*A[J]+C;
- A[J]:=P1 MOD 10;
- C:=P1 DIV 10;
- END;
- END;
-
- PROCEDURE SAVEA;
- BEGIN
- FOR J:=0 TO M DO
- T[J]:=A[J];
- END;
-
- PROCEDURE TERM1;
- BEGIN
- I:=M+M+1;
- A[0]:=4;
- DIVIDEA(I*25);
- WHILE I>3 DO
- BEGIN
- I:=I-2;
- CLEARP;
- P[0]:=4;
- DIVIDE(I);
- SUBA;
- DIVIDEA(25);
- END;
- CLEARP;
- ADJUST;
- SUBA;
- DIVIDEA(5);
- SAVEA;
- END;
- PROCEDURE TERM2;
- BEGIN
- I:=M+M+1;
- A[0]:=1;
- DIVIDEA(I);
- DIVIDEA(239);
- DIVIDEA(239);
- WHILE I>3 DO
- BEGIN
- I:=I-2;
- CLEARP;
- P[0]:=1;
- DIVIDE(I);
- SUBA;
- DIVIDEA(239);
- DIVIDEA(239);
- END;
- CLEARP;
- ADJUST2;
- SUBA;
- DIVIDEA(239);
- SUBT;
- END;
-
- {MAIN PROGRAM}
-
- BEGIN
- CLRSCR;
- WRITELN(' THE COMPUTATION OF PI');
- WRITELN(' -----------------------------');
- 10:WRITELN('INPUT NO. DECIMAL PLACES');
- READLN(M1);
- TI:=TIME;
- M:=M1+4;
- FOR J:=0 TO M DO
- BEGIN
- A[J]:=0;
- T[J]:=0;
- END;
- TERM1;
- TERM2;
- MULT4;
- WRITELN;WRITELN;
- WRITE('PI = 3.');
- FOR J:=1 TO M1 DO
- BEGIN
- WRITE(A[J]);
- IF J MOD 5 =0 THEN WRITE(' ');
- IF J MOD 50=0 THEN WRITE(' ');
- END;
- WRITELN;WRITELN;
- WRITELN('STARTING TIME = ',TI);
- WRITELN('ENDING TIME = ',TIME);
- WRITELN;
- GOTO 10;
- END.